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SOME  IMPLICIT  FINITE  DIFFERENCE  SCHEMES 
FOR  HYPERBOLIC  SYSTEMS 


ABSTRACT 

Implicit  finite  difference  schemes  usually  allow  the  use  of  large  time  steps. 
This  could  be  an  advantage  in  computing  slowly  varying  solutions  to  hyperbolic  sys- 
tems.    The  usual  Crank-Nicholson  scheme  as  applied  to  a  hyperbolic  system  is 
rather  slow.     This  paper  describes  two  modifications  of  the  Crank-Nicholson 
scheme  for  hyperbolic  systems  which  are  implicit  only  in  that  they  require  the  in- 
version of  tridiagonal  matrices.     However,   these  schemes  are  unconditionally 
stable  only  for  positive  definite  systems  (supersonic  flow).     The  results  of  compu- 
tations using  these  schemes  are  described. 
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SOME  IMPLICIT  FINITE  DIFFERENCE  SCHEMES 
FOR  HYPERBOLIC  SYSTEMS 

by 
John  Gary 


1.     Introduction 

In  Section  2  we  will  describe  the  application  of  the  Crank-Nicholson  scheme 
to  hyperbolic  systems  of  equations.     This  method  requires  the  inversion  of  a  block- 
tridiagonal  matrix.     In  order  to  speed  up  the  computation  we  modify  the  scheme  so 
that  we  need  only  invert  scaler  tridiagonal  matrices.     All  of  our  results  are  for 
problems  in  one  space  dimension.     These  methods  can  be  generalized  to  two  dimen- 
sions but  we  have  no  analysis  to  indicate  that  the  generalizations  will  work.     In  Sec- 
tion 3  we  describe  the  numerical  experiments  used  to  check  the  accuracy  and  sta- 
bility of  these  schemes.     The  equations  of  inviscid  hydrodynamic  flow  were  used  in 
these  calculations.     In  Section  4  we  analyze  the  stability  of  these  schemes  using  the 
method  of  von  Neumann.     In  some  cases  we  are  unable  to  explicitly  determine  the 
eigenvalues  of  the  amplification  matrix  so  we  are  forced  to  compute  the  eigenvalues 
numerically. 

2.     Description  of  the  Finite  Difference  Schemes 

In  this  section  we  will  describe  some  finite  difference  schemes  of  an  im- 
plicit nature.     For  purposes  of  comparison  we  will  begin  with  the  usual  Crank- 
Nicholson  scheme  [4].     This  scheme  is  adapted  to  the  nonlmear  problem  by  use  of 
a  "predictor- corrector''  method.     We  denote  the  hyperbolic  system  of  equations  by 
3w/9t  +  A(w)  9w/9x  =  0    where    A(w)    is  a  matrix  with  real,   distinct  eigenvalues. 
If  the  equations  are  those  of  hydrodynamic  flow  we  have 
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where    p,    p,    and  u    denote  density,    pressure  and  velocity  and    y    is  the  ratio  of 

specific  heats.     We  will  use  a  mesh  of  equally  spaced  points    (x  ,  t   )    where 

J     n 

X..,  -  X    =  Ax,     t         -  t     =  At    and    1   <  j   <   M.       We  use  the  notation 
j+1         J  n+1        n  -  -^  - 

w.    =  w(x.,  t   ), 
J  J     n 


w 


w_ 

X 


w(x+Ax,  t)  -  w(x,t) 
Ax 

w(x,  t)  -  w(x-Ax,  t) 
Ax 


w(x+Ax,  t)  -  w(x-Ax,  t) 

W/\      -     — — — 

X  2  Ax 

We  shall  assume  the  boundary  values    w(x     t)    and    w(x     ,  t)    are  known  and 
constant.     Of  course,    this  problem  is  not  properly  posed,    but  where  the  bound- 
ary conditions  are  known  and  constant  this  does  not  seem  to  cause  any  trouble. 

The  first  step  in  the  Crank-Nicholson  method  is  to  predict  the  values 

.       n+1 
of    w.  by  the  use  of  an  explicit  difference  scheme.     We  denote  the  predicted 

values  by    w^,     and  define  them  by 

w''     =    w       -     At  A(w   )  Wy,      . 

X 

n+1 
Then    w  is  defined  by  the  boundary  conditions  and  the  following  equations 


n+1  n 

w  =    w 


n 


^t    ,,w"'  +  w*,,    n+1  n, 

— -  A( )(w^       +  w^) 

2  2  X  X 


2  2 

The  truncation  error  is    0(Ax    +  At   ), 

To  solve  these  equations  we  use  a  method  given  by  S.   Schechter  [5]. 
We  may  write  the  system  as 


-  5  - 


o  .      n+1  n+1  r> .      n+1  ^ 

pA  w.    ,     +    w.         -    /3A.W.    ,     =    D.      , 
J    J+1  J  J    J-1  J 

where    ^  =  At/(4Ax),     A.  =  A((w.       +  wf)/2)    and    D    is  a  vector  defined  by 

«f  J  «l 

D.  =  w.     -    pA.(w.    ,   -  w.    ,)      . 
J         J  J      J+1  j-1 

Thus  the  matrix  of  this  system  is  block  tridiagonal.     The  method  is  analogous  to 

that  used  for  a  scaler  tridiagonal  system.     We  define  matrices    F.    and  vectors 

2         - 1 
G.    by  recursion  as  follows:      F,   =1,     G,   =  D, ,     F.  =  1  +  ^  A.F.   ,  A.    ,,     and 
J  _^  1  1  1         J  J    j-1     j-1 

G.  =  D.  +  ^A.F.    ,D.    ,.      Then  the  solution  is  obtained  by  "backward  substitution" 
J         J  J    J-1    J-1 

n+1      „-l^  n+1      „-l,^         o,      n+1, 

w,^     =F.^G^^,     w.        =  F.      G.  -  ^A.w.^,    . 

M  M     M         J  J        3  J    j+1 

To  use  this  method  we  must  be  certain  that  the  matrices    F.    can  be  in- 

3 
verted  without  difficulty.     If  the  matrix    A(w)    is  constant,     then    A    has  a  complete 

set  of  distinct  eigenvectors.     These  are  also  a  complete  set  for  the  matrices    F.. 

If   a    is  an  eigenvalue  of    A.    and    X.    is  the  corresponding  eigenvalue  of    F.,     then 

2    2  J  J  11       11  2 

1  +^  a   /X.    IS  the  corresponding  eigenvalue  of    F.      .      Therefore    1<  J|F.jj  <  1+a 

for    1   <  j  <  M,     where    a  >  ^j|a|].      If  the  matrices    A    are  all  symmetric,   then 

the  norm  of    F.    satisfies  the  same  inequality.     In  our  case  the  matrices    A.    are 

J  J 

determined  by  equations  (2. 1)  and  are  therefore  neither  symmetric  nor  constant. 

However,   we  had  no  difficulty  in  inverting  the  matrices    F..      An  analysis  by  the 

method  of  von  Neumann  [4]  shows  this  scheme  to  be  unconditionally  stable  and  this 

conclusion  is  supported  by  numerical  experiments  described  in  the  following  section. 

We  will  next  define  a  quasi-Crank-Nicholson  (Q-C-N)  scheme  which  is  not 

centered  in  time.    We  define  the  lower  triangular  matrix    A       to  consist  of  those 

elements  of  A  on  and  below  the  diagonal  with  zero  elements  above  the  diagonal.    We 

define    A       by    A     =  A  -  A    .     Then  the  (Q-C-N)  difference  scheme  is  defined  by 
U  U  L 

w  =    w      -    AtA    (w  )w^        -    AtA    (w  )  w^     .  (2.2) 

To  solve  for  each  dependent  variable  we  need  only  solve  a  system  of  equations 
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whose  matrix  is  tridiagonal.     Thus  the  scheme  is  effectively  explicit.     This  is 
the  reason  for  splitting  the  matrix  into  triangular  parts.     If  the  matrix    A(w) 
is  defined  by  equations  (2.  1);    then  this  scheme  is  unconditionally  stable  for 

supersonic  flow  and  unconditionally  unstable  for  subsonic  flow.     This  will  be 

2 
discussed  in  Section  4.     The  truncation  error  is    0(At  +  Ax   ). 

We  can  also  define  a  quasi-Crank-Nicholson  scheme  which  is  centered 
in  time.     The  first  step  is  to  predict  the  values  of    w  by  the  equations 

,„    ox  *  ri  .     .  ,    n,     n 

(2.3)  w*    =    w       -    AtA(w   )wy.      . 

X 

The  centered  Q-C-N  scheme  is  then  defined  by  the  equations 

,n    A,  n+1  n  At       u      n+1  At       #      #  At      #     n 

(2.4)  w  =    w       -     -  a[w^         -     -  A^w^    -     T  ^    "$     ' 

#  n        # 

where    A     =  A(  (w    +w   )/2).       The  solution  of  this  system  is  obtained  by  solving 

a  tridiagonal  matrix  equation  for  each  dependent  variable.     The  truncation  error 

2  2 

is    0(At    +Ax  ).      In  Section  4  we  will  show  that  this  scheme  is  unconditionally 

stable  for  supersonic  flow  and  conditionally  stable  for  subsonic  flow,    provided 

that  the  matrix    A    is  given  by  equations  (2.  1). 

This  centered  quasi- Crank- Nicholson  scheme  was  applied  to  the  Navier- 
Stokes  equations: 

dp  dp  8u 

9t  9x  dx 

9T  9T  ,        ,v^   3u  7        a^T  4y     ,        ,w3u,2 

r         9x 

2 

9u  au  1   9T  T_  9p     _       4      9_u 

9t  "  9x     ^    7   9x      '*'    7p    9x      ^     3pR   ^    2       ' 

ox 

In  these  equations    R    represents  the  Reynolds  number  and    P      the  Prandtl  num- 
ber.    If  these  equations  are  differenced  by  the  Lax-Wendroff  method,   then  the 
stability  condition  is  approximately 
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At     <     min 


2  -I 
Ax         3Rp  Ax 

16 


[3]. 


There  is  a  parabolic  dependence  on    Ax.      We  have  not  determined  a  stability 
criterion  for  the  centered  Q-C-N  method  applied  to  the  Navier-Stokes  equations. 

However,    enough  computations  were  carried  out  to  show  that  the  stability  con- 

2 
dition  depends  on    Ax      very  weakly,   if  at  all.     Thus  the  stability  condition  is  less 

stringent,   especially  for  values  of    R    which  are  not  large. 

It  might  be  possible  to  use  these  methods  for  problems  in  two  space 
dimensions  by  using  alternating  direction  techniques.     For  example,   suppose  we 
wish  to  solve  the  system 

-—    +    A(w)  ^-    +   B(w)  ^—     =    0      . 
at  9x  ay 

The  finite  difference  scheme  could  be 

#  n  A.  aO     #  ...on  ,    „o     n 

w*    =    w      -    AtA^  wi.    -    AtA^^w^    -    AtB    wyv     , 
L.     X  U    X  y 


w 


n+1 


-    w""    -    AtA^wi    -    AtEN"""^    -    AtB^w#     , 


where    A     =  A(w  )    and    A     =  A(l/2(w    +w  } ).      If  this  scheme  is  stable  at  all,   it 
would  probably  be  stable  only  for  supersonic  flow.     The  fact  that  the  truncation 
error  is  first  order  in    At    may  not  be  a  major  defect,   particularly  if  the  method 
is  used  to  solve  for  a  steady  state  solution.     Perhaps  it  is  possible  to  devise  a 
scheme  which  is  unconditionally  stable  for  both  subsonic  and  supersonic  flow  which 
involves  nothing  worse  than  the  solution  of  scaler  tridiagonal  matrices. 

3.     Results  of  Numerical  Computations 

To  test  these  schemes  we  used  a  hydrodynamic  flow  containing  a  simple 
rarefaction  or  compression  wave.     We  applied  the  difference  schemes  to  equations 
(2.  1).     The  initial  value    u(x,  0)    was  chosen  to  be  constant  in  the  neighborhoods  of 


X      and    X       (the  boundary  points)  and  monotone  in  between.     The  initial  values 
p(x,  0)    and  p(x,  0)    were  computed  from  the  equations  below.     This  produces  a 
simple  wave  moving  on  the  characteristic  with  slope    u  +  c. 

p(x.O)    =pJl.I^i]2/<^-»     , 

o 

p 

p(x,  0)     =   —  p(x,  0^      . 

In  these  equations    p   ,    p  ,    c      represent  the  constant  state  ahead  of  the  wave. 

o       o       o 

The  exact  solution  to  this  problem  is  easily  calculated  [1], 

The  results  of  the  computation  are  shown  in  Table  I.     The  number  of 
points  in  the  first  column  refers  to  the  number  of  mesh  points  used  to  cover  the 
interval  over  which    u(x,  0)    is  not  constant.     This  determines  the  value  of    Ax. 
The  value  of    At    is  given  by    At  "  0.  9  Ax/(  j  uj  +  c)    which  is  the  Courant-Friedrichs- 
Lewy  condition  for  stability.     The  uncentered  Crank-Nicholson  (C-N)  scheme  re° 
fers  to  one  in  which  the  coefficient  matrix  is  not  centered,   that  is: 

,^  ,  ,    n, ,    n+1         n. 
^1  AtA(w  )(w^       +  Wa) 

n+1  n  XX 

.    w — . 

In  one  case  the  quasi-Crank-Nicholson  scheme  was  run  with  an  initial  value 

2 
u(x,  0)    which  had  a  continuous  third  derivative.     In  the  other  cases    u(x,  0)   e   C 

3  3 

but    9   u/9x      was  discontinuous  at  the  head  and  tail  of  the  rarefaction  wave.     If 

2 
the  truncation  error  is    0(  Ax  ),     the  percentage  error  should  drop  by  a  factor  of 

eight  when  the  mesh  spacing  is  halved.     The  percentage  error  given  in  the  table 

is  the  maximum  percentage  error  in  the  variables    p,    p,    u    throughout  the  mesh 

2 
at  100  time  steps.     If  the  truncation  error  is  to  be    0(  Ax  ),     we  must  have 
3 

u  €  C   .      The  flow  is  a  simple  rarefaction  wave  with  Mach  number    M     -  2 

o 

ahead  of  the  wave  and    M     -  1.  75    behind  the  wave.     Other  runs  were  made  with 

M     =0    and    M,   -0.7.      In  this  case  the  error  was  about  half  that  shown  m  the 
o  1 

table.     The  computing  time  shown  at  the  bottom  of  the  table  is  the  time  required 
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Table  I.     Percentage  error  after  100  time  steps  for  the  various  difference  schemes. 


No.   of 
Points 

C-N 
(A  cen,) 

C-N 
(A  not  cen. ) 

Q-C-N 

(cen. ) 

Q-C-N 
(cen. ) 
ueC3 

Q-C-N 

(not  cen. ) 

Lax- 
Wend. 

10 

20 

100 

200 

4.4% 
1.7 
0.055 
0.0087 

4.7% 
1.9 
0.  11 
0.042 

4.1% 
1.6 
0.047 
0.0073 

4.8% 
1.8 
0.039 
0.0050 

8.8% 
5.2 
0.57 
0.15 

2.0% 
0.71 
0.014 
0.0020 

Computation 
time  in  milli- 
seconds per 
mesh  point. 

15 

13 

4.4 

4.4 

1.7 

3.8 

(in  milli-seconds)  to  compute  the  solution  at  a  single  mesh  point.     Results  of  com- 
puting with  the  Lax-Wendroff  scheme  are  shown  for  comparison.     The  equations 
used  in  the  Lax-Wendroff  scheme  are  written  in  conservation  form.     This  accounts 
for  much  of  the  increase  in  accuracy.     In  Table  I,   C-N  denotes  Crank-Nicholson, 
Q-C-N  quasi" Crank- Nicholson  and  cen.   denotes  a  centered  scheme. 

4.     Derivation  of  Stability  Conditions 

We  will  first  derive  a  stability  condition  for  the  centered  quasi-Crank- 
Nicholson  scheme  with  the  matrix    A(w)    given  by  equations  (2. 1)«     We  use  the 
method  of  von  Neumann,   that  is,   we  linearize  the  equations  and  assume  a  pertur- 
bation of  the  form    w     =  w(x.,t  )  =  k     ej 
into  equations  (2.  3)  and  (2,4)  we  obtain 


bation  of  the  form    w"  =  w(x.,t  )  =  k*^  exp  (iwx.).      Substituting  this  value  for    w 

J     n  J  J 


n+1 


I    - 


t^T-J 

-i 

I  - 

'^  a1 

2    ^U_ 

i^A  j   k'^ 


Rk 


10  - 


where    /3  =  (At  sinujAx)/Ax.      If  we  let    a  =  ^u/2,     then  the  amplification  matrix 
R    is  given  by    R  =  I  +  2q'F/(1  -\-ia)    where  the  matrix    F    is  given  by 


F    = 


—  1 


a 


u 


-(i  + 


a 


M 


2        2 
i{a    /M    -  1) 

pud  +  ior) 


ip(l  - 

ia) 

u 

ipu(l  - 

la) 

M 
a      (1   -  io-) 


M 


2   (1  +  ior) 


-  1 


Here    M    denotes  the  Mach  number.     If  we  denote  the  eigenvalues  of    F    by   M• 
(j  =  1,  2,  3),     then  the  eigenvalues  of    R    are    (1  +  ior  +  20?^  ) /d  +  i«)-      L-et 
z.  =  1  +  iff  +  2Qf/Lt.;      then  by  expanding  the  characteristic  equation  of    F    to  find 
an  equation  for   /i .    we  obtam    z     =  1  -  ior,     (z.  +  Q  )(z.  +  Q  )  =  -  Q      (j  =  1,  2) 


J 


where 


J 


J 


Q. 


=    Qfl     + 


2a 
M 


-     1 


Q. 


ai    + 


2a^{l 


iff) 


M  (1  +  iff ) 


-     1 


Q. 


4ff^(l  -  ff^/M^)(l  -  Iff) 
M^(l  +iff) 


-1    -   "^  ^^1  +  Q.,)  =  (1  -  i'*)(2  -  4ff    /(M  (1  +ff    ) )    and    z    z     =  (1 


Then     z,  +  z„  =  -  (Q. 


-  iff)' 


Let    z.  =  a.d  -  iff),     then    a  a„  =  1    and    a    +  a„  =  2  -  4ff   /(M  (1  +ff  ) ).      The 
J         J  12  12 

eigenvalues  of    R    are    r     =  (1  -  iff)/(l +iff)    and    r.  =  z./(l+iff)    (j  =  1,  2).      Thus 

•^  J        J 

|r    1=1,     and    |r.j-|a.|     (j  =  1,  2).      Therefore    |r.|<l    (j  =  1,  2,  3)    if  and 

only  if   ff     <   M  (1+  ff   ).      This  follows  from  the  product  and  sum  relations  for 

a..      Therefore  we  have  unconditional  stability  if    M  >   1,      If    M  <   1,     then  the 

222  22  22 

inequality   a     <   M  (1+ff   )    reduces  to    ^    c     <   4  +3   u      or    At  <   2 Ax 


/V7^ 


where    c    is  the  velocity  of  sound.     This  completes  the  stability  analysis. 
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We  will  now  analyze  the  stability  of  the  uncentered  scheme  defined  by  equa- 
tions (2.  2).     Again  we  use  the  von  Neumann  method.     The  amplification  matrix    R 
is    R  =  C'^C      where    C     =  I  +  ^A    ,     C^  I  -  0A    ,     and   0  ■- i  At  sin  u  Ax/Ax.      We 
first  assume  that    A(w)    is  constant,    symmetric  and  positive  definite.     We  can  then 
prove  that     ||r||    <   1.      The  proof  is  almost  exactly  the  same  as  that  given  by 
Ostrowski  to  prove  that  the  Gauss-Sidel  iterations  converge  [2].     Assume,   for 
simplicity,   that  the  order  of    A    is  three.     Let    e     =(1,0,0),     e     =  (0,  1,  0)    and 
e     =(0,0,  1).      Let    x      be  an  arbitrary  three-dimensional  vector,   and  let    [x] 
denote  the  i-th  component  of  a  vector    x.      Let    x     =  x    +  o  e  ,     x     -  x    +  "n^o    ^^'^ 
X     =  X     +  car   e      where    a     is  chosen  such  that    [C  x.].  =  [C  x.      ].      (j  =  l,  2,  3).      It 

is  easy  to  show  that    C,x„  -  C„x   .      We  continue  this  process  to  obtain  a  sequence 
■^  1   3         2  o 

of  vectors    x.    such  that    C  x      =C  x         .       We  will  show  that    jx.|  -0    (for  0^0) 
J  l_j'J  2    oj-o  J 

which  can  only  occur  if    j|  C~    C    ||    <   1.      Thus  we  will  have  proved  stability.     Since 
A    is  symmetric  and  positive  definite  we  can  define  a  norm  for  complex  vectors    x 
by    llxll   =x*Ax.      We  have    x..,   =x.+a.      e,     where    1  <  k  <   3.      A  little  algebra 
will  yield 


llxjil     -     l|x..J|     =    2Rei«.[Ax..,]J 


+    or.or.a,  , 

3    J   kk 


From  the  definition  of   x.    we  have    [C  x       -  C  x.]     =  0. 


From  this  we  have 

[i3Ax.],    =a.(l-^a,,).      Therefore    ||  x       ||   -  ||  x  ||   =  -a  or  a       <  0.      This  in  turn 
3  k        J  kK  jTi  J  J    J    t^*^ 


yields  the  convergence  of    ||x.      ||     which  shows  that    a.  —  0.      Since 

r3Ax  1     =  «.(1  -^a,  ,  )    we  can  prove  that    x.  —  0. 
j  k        J  kk  J 

If  the  matrix    A    is  given  by  equations  (2.  1)  we  are  no  longer  able  to  obtain 
a  stability  condition  by  analytic  means.     However  we  can  compute  the  eigenvalues  of 
the  amplification  matrix  numerically  and  thus  determine  a  stability  condition.     The 
difference  equation  is 

n+1  n  .     .     /    n.     n+1  .^  .     ,    n      n 

w  =    w      -    At  A^  (w  )  w/v        -    At  A_,(w  )  wa     . 

L  X  U  X 

If   A(w)    is  defined  by  equations  (2.1),  then  the  amplification  matrix  for  this  scheme 
is 
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Figure  1.     Stability  of  the  Q-C-N  scheme  (uncentered) 
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R    = 


1 
1  +  ia 

0 


0 
1 


1  +  iof 

-  iff 
pd+ior)' 


i/3p 
1  +  la 

jffTP 
1  +  ice 


-  ff  7p 
Pd+iof)' 


1  +  iof 


where    /3  =  (At  sin  u  Ax) /Ax    and   or  -  /3u.      The  eigenvalues  of    R    are  given  by 
H  ./(I  +ia)    where   ^„  =  1    and   p  .    (j  =  1,  2)    are  the  roots  of  the  following 
quadratic 


-     (2-6)^    +    1     =    0 


6    = 


r2   2 
P   c 

1  +  ior 


Note  that  the  eigenvalues  of    R    depend  only  on    a  =  /3c    and    M  =  u/c.      We  let 
X    denote  the  absolute  value  of  the  largest  eigenvalue  of    R.      In  Figures  1  and  2 
X    is  plotted  against    a    for  various  values  of  the  Mach  number    M.      These 
graphs  indicate  that  the  scheme  is  unconditionally  stable  for    M  >   1    and  uncon- 
ditionally unstable  for    M  <   1.      Numerical  computations  were  performed  using 
this  scheme  as  described  in  Section  3.     The  flow  was  a  rarefaction  wave  with 
M  =  2.  0  ahead  and    M  =  1.  75    behind  the  wave.     If  the  value  of    At    was  made 
equal  to  ten  times  the  value  of  Courant-Friedrichs-Lewy  (i.  e.  , 
At  =  10  Ax/(juj  +  c) ),     then  there  was  no  sign  of  instability  out  to  44  time  steps. 
Computations  were  also  performed  with  a  subsonic  rarefaction  wave  with 
M  =  0,0    ahead  and    M  =■  0.  7    behind  the  wave.     With    At  =  2  Ax/(  |uj  + c)    the 
density  became  negative  at  61  time  steps  and  when    At  =  Ax/(juj  +  c)    the  density 
was  negative  at  377  time  steps.     For  small  values  of    At    the  largest  eigenvalue 
of  the  amplification  matrix    R    is  very  close  to  the  unit  circle,   therefore  in- 
stability is  slow  to  develop. 

The  eigenvalues  of    M    drop  off  quite  rapidly  for  supersonic  flow  accord- 
ing to  Figure  2.     This  suggests  that  the  scheme  might  produce  reasonable  results 
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even  if  the  flow  contains  a  shock.     We  tested  this  scheme  on  a  flow  which  was  ini- 
tially an  isentropic  compression  wave.     Eventually  the  compression  wave  will 
produce  a  shock.     This  provides  a  test  of  how  the  difference  scheme  will  react 
to  the  shock.     We  ran  the  Lax-Wendroff  scheme  on  the  same  flow  for  a  compari- 
son.    The  initial  compression  wave  had  a  Mach  number  of  2.  53  behind  and  1.  5 
ahead,   which  produced  a  very  strong  shock.     The  computed  values  of  pressure  are 
plotted  in  Figure  3.     Both  schemes  were  run  to  300  time  steps  (at    Ax  =  0.  05).      By 
printing  out  the  computed  value  of  pressure  every  40  time  steps  from  100  to  300 
time  steps  and  observing  where  the  shock  is  located  we  can  compute  the  shock 
velocity.     For  both  schemes  the  shock  velocity  is  constant  within  the  accuracy  of 
measurement.     We  know  the  state  ahead  of  the  shock;    thus  we  can  compute  the 
state  behind  the  shock  from  the  shock  velocity  and  the  Rankine-Hugoniot  relations. 
The  results  are  given  in  the  table  below. 


Lax-Wendroff 

Quasi-C.  -N. 

Observed 

Computed 
from  R-H 

Observed 

Computed 
from  R-H 

Shock  Mach  number 
Pressure  ratio 
Density  ratio 
Velocity  behind 

3.00 
10.2 
3.84 
3.70 

10.3 
3.85 
3.72 

2.33 
10.7 
5.4 
3.65 

6,17 
3.12 
3.11 

When    Ax    is  reduced  to  0.02  the  results  from  the  quasi-C.N.   scheme  do 
not  improve.     The  observed  shock  speed  is  2.  34  for    Ax  =  0.02.     The  computed 
value  of  pressure  for    Ax  =  0.  02  is  plotted  in  Figure  3.     This  lack  of  improvement 
is  to  be  expected  since  the  state  immediately  behind  the  shock  is  nearly  constant. 
Thus  the  results  from  the  Q.  -C. -N.   scheme  are  apparently  worthless,   although 
the  curves  in  Figure  3  appear  reasonable  upon  casual  inspection. 
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Figure  3.     Pressure    p    vs.   distance    x    at  150  time-steps 
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Q-C-N  (uncentered) 
(Ax  =  0.  02,    displaced 

vertically  one  unit) 


Figure  3    (Continued) 
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